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ABSTRACT 

We compute statistical properties of weak gravitational lensing by large-scale structure in 
three Cold Dark Matter (CDM) models: two flat models with (fi ,A ) = (1,0) and (0.3,0.7) 
and one open model with (f2o,Ao) = (0.3,0), where fio and Ao are the density parameter and 
cosmological constant, respectively. We use a Particle-Particle/Particle-Mesh (P 3 M) iV-body 
code to simulate the formation and evolution of large-scale structure in the universe. We 
perform 1.1 x 10 ray-tracing experiments for each model, by computing the Jacobian matrix 
along random lines of sight, using the multiple lens-plane algorithm. From the results of these 
experiments, we calculate the probability distribution functions of the convergences, shears, and 
magnifications, and their root-mean-square (rms) values. We find that the rms values of the 
convergence and shear agree with the predictions of a nonlinear analytical model. We also find 
that the probability distribution functions of the magnifications /x have a peak at values slightly 
smaller than [i = 1, and are strongly skewed toward large magnifications. In particular, for 
the high-density (CIq = 1) model, a power-law tail appears in the distribution function at large 
magnifications for sources at redshifts z s > 2. The rms values of the magnifications essentially 
agree with the nonlinear analytical predictions for sources at low redshift, but exceed these 
predictions for high redshift sources, once the power-law tail appears. 

We study the effect of magnification bias on the luminosity functions of high-redshift 
quasars, using the calculated probability distribution functions of the magnifications. We 
show that the magnification bias is moderate in the absence of the power-law tail in the 
magnification distribution, but depends strongly on the value of the density parameter £Iq. In 
presence of the power-law tail, the bias becomes considerable, especially at the bright end of the 
luminosity functions where its logarithmic slope steepens. We present a specific example which 
demonstrates that the bias flattens the bright side logarithmic slope of a double power-law 
luminosity function. 



Subject headings: cosmology: observations — cosmology: theory — gravitational lensing 
large-scale structure of universe — quasars: general 



1. INTRODUCTION 

It is well known that the apparent brightness of distant sources are gravitationally affected by 
the inhomogeneous distribution of matter in the universe, an effect called the lensing magnification. 
Consequently, the apparent brightness and intrinsic luminosity of a distant source are no longer simply 
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related by the luminosity distance-redshift relation in a smooth Friedmann universe (Weinberg 1972; 
Schneider, Ehlers, & Falco 1992, hereafter SEF). The relation depends on how light rays are being lenscd 
by density inhomogcneitics in the universe, as they propagate from the source to the observer. Since the 
pioneering works of Kristian & Sachs (1966) and Gunn (1967), numerous studies of the effect of lensing 
magnification in an inhomogeneous universe have been published. Two different approaches have been used 
in these studies to compute the statistical properties of lensing magnifications in inhomogeneous universes 
with realistic matter distributions: the power spectrum approach, which is analytical, and the numerical 
approach, which combines TV-body simulations of structure formation with ray-tracing simulations. 

The power spectrum approach was first developed by Gunn (1967), who showed that the root-mean- 
square (rms) fluctuation in the apparent brightness of distant sources (or cquivalcntly in the lensing 
magnification) can be expressed as a radial integral of the density autocorrelation function. Gunn's 
original approach was eventually modified by Babul & Lee (1991) to reflect substantial advances in our 
understanding of large-scale structure formation. These authors gave a quantitative estimate of the rms 
fluctuations for modern cosmological models such as the Cold Dark Matter (CDM) model. Fricman (1997) 
rederived the rms fluctuations using the Limber's equation in Fourier space (Kaiser 1992, 1998), and found 
that the dispersion in the lensing magnification can be larger than 0.1 for sources at redshift z s = 1, but 
depends strongly on the background cosmological model. The skewness of the probability distribution of 
the lensing magnifications was first calculated by Nakamura (1997) using the quasi-linear theory of density 
fluctuation. 

The first study based on the numerical approach was done by Jaroszyfiski et al. (1990), who used 
a Particle-Mesh (PM) iV-body code to simulate the formation and evolution of large-scale structure in a 
CDM universe. The multiple lens-plane algorithm developed by Blandford & Narayan (1986) and Kovner 
(1987) was then used to follow the evolution of ray bundles propagating through the inhomogeneous matter 
distribution in the simulated universes. The force resolution of the PM simulations was l/i _1 Mpc, in 
present units (where h is the Hubble constant H in units of 100kms _1 Mpc _1 ), not sufficient to include 
lensing effects by small-scale nonlinear structures like galaxy clusters. They found that the rms values of the 
lensing magnifications by large-scale structure are very small even for sources at redshift z s — 5. Since then, 
this method has been improved by several authors (Bartelmann & Schneider 1991, 1992; Jaroszyhski 1991, 
1992; Wambsganss et al. 1995, 1997; Wambsganss, Cen & Ostrikcr 1998; Premadi, Martel & Matzner 1998; 
Tomita 1998a, 1998b). Thanks to their effort and recent rapid developments of computational techniques as 
well as computing power, the length resolution of these algorithms has significantly improved. For instance, 
Wambsganss et al. (1998) used large PM TV-body simulations combined with a convolution method and 
have achieved an effective resolution of 10/i _1 kpc. 

Jaroszyfiski (1991, 1992) extended the length resolution of the method to subgalactic scales by 
combining numerical simulations of large-scale structure formation with a Monte-Carlo method for locating 
galaxies inside the computational volume. Galaxies were modeled as isothermal spheres with velocity 
dispersions and core radii determined by empirical relations depending on morphological types. Ray-tracing 
simulations were then performed, taking into account the combined effects of the large-scale structure 
and the galaxies. Later, Premadi et al. (1998) improved this method significantly, first by replacing the 
Zcl'dovich algorithm used by Jaroszyfiski for simulating large-scale structure formation by a more accurate 
Particle-Particle/Particlc-Mcsh (P 3 M) iV-body code, and then by taking into account, when assigning 
morphological types to galaxies, the observed morphology-density relation (Martel, Premadi, & Matzner 
1998 and references therein). 

In this paper, we use the multiple lens-plane algorithm to study the effect of lensing by large-scale 
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structure on quasar luminosity functions. We focus on CDM models with a COBE-normalized tilted power 
spectrum (Bunn & White 1997). The P 3 M simulations for these models were provided by the Texas 
P 3 M Database (Martel & Matzncr 1999). Each model is characterized by the value of 4 parameters: the 
density parameter Q 07 cosmological constant A , Hubble constant H , and rms density fluctuation er 8 at 
scale 8/t Mpc (the tilt n of the power spectrum is a dependent parameter). We consider three particular 
models, an Einstein de Sitter model (E-dS) with Q — 1 and A = 0, an open model (O) with fl = 0.3 and 
A = 0, and a flat model (A) with fl = 0.3 and Ao = 0.7. For each model, we follow the propagation of 
light ray bundles through the matter distribution. We are primarily interested in the statistical properties 
of lensing by the matter inhomogeneities in the universe, and not in rare events such as multiple imaging of 
distant quasars. Consequently, we focus on the lensing effects by the large-scale (> 0.1/i _1 Mpc) structures 
and make no attempt to study lensing effects caused by individual galaxies. 

In order to obtain good statistics, we perform a total of 1.1 x 10 7 ray-tracing experiments for each 
model. Each experiment consists of computing the Jacobian matrix along a random line of sight. Having 
such a large number of experiments allows us to study in detail not only rms values, but also the probability 
distribution functions of the lensing properties, i.e., the convergence, shear, and magnification. We compare 
these rms values with predictions based on the power spectrum approach. This enables us to test the 
validity of the power spectrum approach and its limitation. 

Using the probability distribution functions of the lensing magnifications obtained from the experiments, 
we study the effects of the magnification bias on quasar luminosity functions. This problem was first 
studied by Sanitt (1971), and later by many different authors (e.g., Turner 1980; Avni 1981; Canizares 1982; 
Peacock 1982; Victri 1985; Ostrikcr & Victri 1986; Schneider 1987b, 1987c, 1992; see also SEF, chapter 12). 
A crucial difference between these previous studies and ours is that we compute the probability distribution 
in universes with realistic matter distributions, that originate from the growth of primordial density 
fluctuations with a specific power spectrum and normalization, and are simulated using a state-of-the-art 
Af-body code. In previous studies, the probability distribution was calculated by considering simple matter 
distributions such as randomly distributed point masses and/or isothermal spheres. These studies revealed 
that the magnification bias can have a considerable effect on the luminosity functions, especially for bright 
quasars (Ostriker & Vietri 1986; Schneider 1992). However, they cannot be used to make predictions 
about the effect on the magnification bias in specific cosmological models, since this requires a realistic 
representation of the large-scale structure in the universe. 

We demonstrate the importance of the bias on quasar luminosity functions and show its dependence 
on the mean matter density of the universe. However, we do not attempt to derive a precise quantitative 
estimate of the bias on the observed luminosity functions. This would require a precise knowledge of both 
the observed luminosity functions and the probability distribution of the lensing magnifications. However, 
uncertainties in the quasar luminosity functions are still considerable (e.g., Boyle, Shanks, & Peterson 1988; 
Hartwick & Schade 1990; Warren, Hewett, & Osmer 1994; Hawkins & Veron 1995; La Franca & Cristiani 
1997). Furthermore, the magnification distributions we calculate are those for a point like source, the effect 
of finite source size on the magnification probability becomes important at large magnifications (Schneider 
1987b; Schneider & Weiss 1988a, 1988b). Because of these limitations, we decided to approximate the 
intrinsic quasar luminosity functions using either a single or a double power-law model. Although it is 
not clear whether the power-law models provide an accurate description of a real luminosity function, we 
believe that these models are sufficiently realistic (especially the double power-law model) for a qualitative 
study of the essential properties of the magnification bias. 

The remainder of this paper is organized as follows: In §2, we describe the cosmological models, and the 
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numerical methods used for simulating both the large-scale structure formation and the light propagation. 
In §3, we present the results of the experiments and compared them with predictions based on the power 
spectrum approach. In §4, we discuss the effects of magnification bias on the quasar luminosity functions. 
Summary and conclusion are presented in §5. In Appendix A, we summarize the basic equations used in 
the power spectrum approach. 



2. DESCRIPTIONS OF THE COSMOLOGICAL MODELS AND NUMERICAL 

METHODS 

2.1. Cold Dark Matter Models 

We consider three tilted CDM cosmological models. Each model is characterized by its values of Oo, 
Ao, Hq, and ag, with the primordial exponent n being a dependent parameter (Martel & Matzner 1999; 
Premadi et al. 1999). Table 1 lists the models and gives their parameters. The power spectrum is given by 

P(k) = 2^(^X +n 5lk n THk), (1) 



H, 



o 

where c is the speed of light. The transfer function T(k) for CDM models is given by Bardeen et al. (1986) 
as follows, 

T(q) = ^^^ [l + 3-89? + (16.1q) 2 + (5.46g) 3 + (6.71g) 4 ]- 1/4 , (2) 
where q is defined by 

« = (mp^i)"" 172 ^^ 2 )" 10 -^ ( 3 ) 

a = a -Wfio a -(Wfio) 3 7 ( 4 ) 

ai = (46.9r> / l 2 ) a670 [1 + (32.1f!o/i 2 )- - 532 ] , (5) 

a 2 = (12.0f>ofe 2 ) a424 [1 + (45.0n ^ 2 r - 582 ] , (6) 

(Hu & Sugiyama 1996), where 62.7 is the cosmic microwave background temperature in units of 2.7K, and 
£!bo is the contribution of baryons to the density parameter. In all models, 62.7 and f2eo were set equal 
to 1 and 0.015ft- -2 respectively. The density perturbation 8r at horizon crossing is obtained by fitting the 



Table 1. Summary of model parameters 



Model fl A flBoh 2 h <r$ n 



E-dS 1.0 0.0 0.015 0.65 1.2 0.8506 
O 0.3 0.0 0.015 0.75 0.85 1.1748 
A 0.3 0.7 0.015 0.75 0.90 0.8796 
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COBE 4-year data (Bunn & White 1997), as follows, 



10 5, = , l.95fi -°-" 191nO ^ 017fl e -("+ - 14 » 2 ) , A = 0; 

H M.94r» - - 785 - - 051nno e -(°- 95 "+ - 169 " 2 ), A = l-fi ; 



where n = n — 1 . 



2.2. The P 3 M Algorithm 

The simulations of large-scale structure formation were provided by the Texas P 3 M Database (Martel 
& Matzner 1999). All simulations were performed using a P 3 M TV-body code (Hockney & Eastwood 1988) 
with 64 3 particles, in a computational cubic box with triply periodic boundary conditions. The forces on 
particles are computed by solving Poisson's equation on a 128 3 cubic lattice using a Fast Fourier Transform 
method. The forces at short distance are corrected by direct summation over pairs of particles separated 
by less than some cutoff distance r e equal to a few grid spacings. The algorithm reproduces accurately the 
Newtonian interaction between particle down to the softening length rj, which is set equal to a fraction of 
the grid spacing. In all simulations, the comoving lengths of the size Lbox of the computational box and the 
softening length rj were equal to 128 Mpc and 300 kpc, respectively. Hence, the algorithm has a dynamical 
range of 427 in length. Notice that in this paper, we express distances in units of /i -1 Mpc. In these units, 
the values of Lbox and rj vary among models. These values, and also the mass per particle, are listed in 
Table 2. Three simulations were performed for each cosmological model. For each model, the simulations 
differ only in the choice of random phases in the initial conditions. All simulations start at an initial redshift 
of z = 24, and end at z = 0. 

We calculated the two-point correlation function of the particles at redshifts z = 0, z = 1, and z = 2 
using a direct estimator (Hockney & Eastwood 1988), 

« r > " SB&v " < 8 > 

where N p is the number of particles inside a spherical shell of inner radius r — Ar/2 and outer radius 
r + Ar/2 centered around a central particle, SV is the volume of this shell, N c is the number of particles 
taken as centers, and n is the mean number density of the particles. The results are plotted in Figure [j]. 
For comparison, we plot analytical estimates of the two-point correlation function, based both on linear 
perturbation theory and on the nonlinear power spectrum fitting formula of Peacock & Dodds (1996), and 



Table 2. Parameters in P 3 M simulations 



Model 


Box Size, ibox 


Particle mass 


Softening length, rj 




(ft- x Mpc) 


ih- l Mv) 


(/i- x Mpc) 


E-dS 


83.2 


3.608 x 10 11 


0.195 


and A 


96.0 


1.249 x 10 11 


0.225 
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Fig. 1. — Two-point correlation functions 
compared with those predicted by the linear 
theory and (a) is for E-dS model, (b) is for O 
model and (c) is for A model. The filled symbols 
represent the two-point correlation functions 
calculated from the particle positions in the 
P 3 M simulations. Filled circles, squares, and 
triangles corresponds to redshifts z = 0, z = 1, 
and z = 2, respectively. Solid lines show the 
correlation function derived from the non-linear 
power spectrum, using the fitting formula by 
Peacock & Dodds (1996). Dashed lines show the 
correlation functions predicted by linear theory. 
The redshifts are z = 0, 1 and 2 from top to 
bottom, respectively. The dotted line shows the 
observed galaxy two-point correlation function. 
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also the observed galaxy 2-point correlation function at the present epoch, 




) 



-1.77 



r 



(9) 



5.4/i- x Mpc 



(Peebles 1993). As we see in Figure [j] the slope the correlation functions drop significantly at short range, 
for r < rj. Martel (1991) pointed out that this flattening of £ is due to the softening of the force, which 
reduces the number of particle pairs with separations r <J rj, and increases the number of pairs with 
separations r > rj. For the E-dS model, the correlation functions obtained from the P 3 M simulation are 
in good agreement with the nonlinear predictions. For the O and A models, the nonlinear predictions 
are slightly smaller than the correlation functions obtained from the simulations at intermediate pair 
separations (rj < r < lft, _1 Mpc). In §2.3, we examine the statistical properties of the gravitational lensing 
by these simulated matter distributions, and compare the rms values of the lensing convergences, shears, 
and magnifications with the predictions of the power spectrum approach using Peacock & Dodds' fitting 
formula for the nonlinear power spectrum. Therefore, the reader should keep in mind that, for the O and A 
models, Peacock & Dodds' fitting formula slightly underestimates the correlation function at short range, 
?' < IhT Mpc, and thus we can expect that the rms values of the convergences, shears, and magnifications 
predicted by the power spectrum approach will be also underestimated. 



We use the multiple-lens plane algorithm to follow the evolution of light rays traveling through the 
simulated matter distributions. The P 3 M simulations provided particle distributions in cubic boxes at 
various redshifts, from the initial redshift z = 24 to to present. By combining these boxes, we can represent 
the matter distribution in the universe along a line of sight extending from the observer to the source.^] For 
sources at redshift z s = 3, this requires 36 boxes for the E-dS model, 38 for the O model, and 47 for the A 
model. In order to eliminate spurious correlations between the large-scale structure in adjacent boxes, we 
combine boxes from different simulations, so that cubic boxes which are directly joined to each other are 
chosen from different simulations. 

In the standard multiple-lens plane algorithm, the matter content of each box is projected onto a single 
plane perpendicular to the line of sight. We increase the accuracy of the algorithm by dividing the cubic 
box into four rectangular subboxes of size Lb ox x ^box x (Zbox/4), and projecting the matter content of 
each subbox onto a plane, thus increasing the number of lens planes by a factor of 4. 

The deflection potential ip l on the i-th lens plane is related to the surface mass density fluctuation 
5T.i(x) — Ej(a?) — (Ei) on that plane by 



where x — (x\,X2) is the position vector in the plane, and xq is an arbitrary cutoff length (its value is 
irrelevant, since only derivatives of (j) 1 have a physical signification). This equation can be rewritten in the 



3 This approach ignores the possible existence of structures larger than Lbox = 128 Mpc. This does not effect on our 
experiments in any significant way, because in the CDM universe, fluctuations at that scale have a very small amplitude even 
at the present time, and their lensing effect is negligibly small. 



2.3. The Ray- Tracing Experiments 




(10) 
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form of a two-dimensional Poisson equation, 

V 2 ip i = 2SE i . (11) 

We solve this equation numerically on each lens plane by first computing the surface density on a 
512 x 512 square lattice from the particle positions, using the Triangular Shaped Cloud (TSC) assignment 
scheme (Hockney & Eastwood 1988, §5.3), and then inverting equation ([□]) using a Fast Fourier Transform 
method (see, e.g., Premadi et al. 1998). Notice that the grid spacings are 0.1625ft.~ 1 Mpc for the E-dS model 
and 0.1875/i _1 Mpc for the O and A models, and are slightly smaller than the softening length -q of the P 3 M 
calculations. The evolution equation of the Jacobian matrix in the multiple- lens plane algorithm is given by 

where Dj is the angular diameter distance between the observer and the j-th lens plane, is the angular 
diameter distance between the i-th, j-th. lens plane, and is an optical tidal matrix defined by, 

n = ( til tl2 \ = f 5 ^ + Utll ~ t22) tl2 

1 \tl2 t22 ) \ tl2 1(^-^22) 

with commas denoting differentiation with respect to x. The last equality was obtained using equation ([0]). 
We compute the second derivatives of ip l on grid points using a standard finite difference formula. 

In general, equation ([l2]) is not an explicit equation for since the summation contains the tidal 
matrix , which must be evaluated along the light ray path. One first has to solve the multiple gravitational 
lens equation in order to compute the location of the light ray on each lens plane and then evaluate Ui at 
that location. However, the deflections of light rays are very small; the deflection angle caused by structures 
larger than 0.1/i _1 Mpc is at most a ~ 1CP 5 . The distance between the actual and unperturbed positions of 
a light ray will be much smaller than the grid spacing of the lattice used for calculating the two-dimensional 
potential and its derivatives. To give an actual example, the transverse distance from the unperturbed ray 
position is at most l± ~ L\> OK a < 10 _3 /i _1 Mpc for a light ray propagating a comoving distance of Lbox- 
We can therefore treat light rays as if they were moving along straight lines. Q This greatly simplifies the 
algorithm. Not only the computation of light ray trajectories becomes unnecessary, but in addition we can 
choose lines of sights that go through the grid points, where the surface mass density fluctuation 5Y<i and 
the derivatives ip\i, 1/^22 an( i V^ia °f the potential are known, thus eliminating the need to interpolate U,; 
from the grid points to the location of the rays. Since the lattice is composed of 512 x 512 = 262 144 grid 
points, there are 262 144 possible lines of sight to choose from. Each experiments consists of selecting one 
line of sight, and solving equation (|lj) recursively along that line of sight for all values of A.,-, starting with 
Ai = I (the identity matrix). It may seem that the number of possible experiments is limited to 262 144, 
but actually, since the P 3 M simulations use periodic boundary conditions, each cubic box can be given 
a random shift, or, equivalently, any grid point in a cubic box can be chosen, independently of the ones 
chosen in other cubic boxes, as long as the same grid point is chosen on the four lens planes within each 




4 We could not make this simplification if we were interested in the deformation and/or multiple imaging of sources with a 
finite size. In that case, it is the relative difference between the deflection angles of nearby rays that matters, and those can be 
quite large. 
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box. The total number of possible experiments is therefore 262 144 box , where A^ox is the number of cubic 
boxes along the line of sight. 

We decompose the Jacobian matrix as 

A= ( 1-^-71 -7 2 — ) f (14) 



-72 + w 1 — 5k + 71 

where 5k represents a fluctuation of convergence from that in smooth Friedmann universe, I7I = (7I +7|) 1 ^ 2 
is the amplitude of shear of a light ray bundle, and uj is the rotation angle of the beam. The image 
magnification factor of a point like source is given by the inverse of the determinant of the Jacobian matrix, 

/i= |detA| = |(1-<5k) 2 -7 2 -w 2 |' ^ 



3. RESULTS 

3.1. Lensing Convergence and Shear 

We first study the statistical properties of the lensing convergence 5k and shear 7. For this purpose we 
perform 10 6 experiments for each model. The probability distributions of 5k and 7 are shown in Figure || 
and|^, respectively, for sources at redshifts z s = 1, 2, and 3. The probability distributions of both 5k and 7 
broaden when the source redshift increases, and those in the E-dS model (top panel) are considerably broad 
compared with those in the O (middle panel) and A (bottom panel) models. The probability distributions 
are slightly broader in the A model than in the O model. These properties arc consistent with analytical 
predictions, as we will demonstrate below. The probability distributions of 5k have a peak at some negative 
value, and are skewed toward positive values. This behavior is prominent in the E-dS model. These 
distributions are a reflection of the mass distribution in the universe: in most of the universe, the matter 
density is smaller than the average value, while most of the matter is concentrated into compact structures 
such as galaxies and clusters of galaxies. Thus, for the majority of the random lines of sight, light ray 
bundles converge less compared to those propagating through a smooth Friedmann universe. The average 
values of the fluctuations of the lensing convergence in our simulations are found to be consistent with zero, 
i.e., \(5k}\ < 10~ 4 . There is a minimum value of the fluctuation of the lensing convergence, which occurs 
when the beam propagates through empty space. It is given by 

x , D DR (z) 

dKmin = 1 - , , , (16) 
D F (z) 

(SEF, chapter 11; Hamana 1998), where Ddr{z) is the so-called Dyer-Roeder distance (Dyer & Roeder 
1972, 1973) with the dumpiness parameter a — (Fukugita et al. 1992), and Dp(z) is the standard angular 
diameter distance. The values of <5K m ; n are given in Table 3 for the three cosmological models considered 
in this paper. Comparing the lower cutoff of the probability distributions found in Figure |^ with these 
minimum values, we find that a value of 5k near 5K m i n is unlikely to occur in our ray-tracing experiments, 
especially for sources at a high redshift. This is consistent with the results of Futamase & Sasaki (1989), 
which are based on a simple analytical approach. 
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Fig. 2. — Probability distributions of fluctuations 
of lensing convergence, 5k. The solid lines, short- 
dashed lines, and long-dashed lines correspond to 
sources at redshifts z a = 1, 2, and 3, respectively, 
(a) E-dS model; (b) O model; (c) A model. 



Fig. 3. — Probability distributions of the lensing 
shear of images. Lines have the same meaning as in 
Figure |. (a) E-dS model; (b) O model; (c) A model. 
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In Figure ||, we plot the rms values a K of the convergence and <r 7 of the shear, defined by 



A T , 



cxp 



— E 



1/2 



(17) 
(18) 



where -/V oxp is the number of experiments. Notice that cr K and ct 7 are expected to be equal, based on 
analytical arguments (Jain & Seljak 1997). Instead, we find that the rms values of the shear computed 
from our ray-tracing experiments are slightly but systematically smaller than those of the convergence. 
However, this difference is probably insignificant, because the shear computed from the experiments might 
be underestimated, whereas the convergence is probably accurate. The shear is mainly caused by the 
traceless symmetric part (Weyl part) of the optical tidal matrix JJ i: which depends upon second derivatives 
of the two-dimensional potential ip 1 . These second derivatives can be underestimated somewhat because of 
the finite resolution of the lattice used for solving equation (|l"l|). The lensing convergence is mainly caused 
by the surface mass density fluctuation <5£j, which is directly computed from the particle positions. 

In each panel of Figure ||, we also plotted three curves, which represent various analytical estimates of 
a K and <r 7 (see Appendix A); the solid curves represent predictions based on the power spectrum approach 
(Kaiser 1992; Bernardeau, van Waerbeke & Mellier 1997; Jain & Seljak 1997; Nakamura 1997), using the 
nonlinear power spectrum by Peacock & Dodds (1996) (hereafter, we refer to this as the nonlinear method), 
dashed curves represent the predictions obtained from a top-hat filtered density field with a smoothing 
scale of R s = 0.1/i _1 Mpc (the filtered nonlinear method), comparable to the softening length of the P 3 M 
simulations, and long-dashed curves represent the predictions of linear perturbation theory (the linear 
method). 

The nonlinear method predicts that the rms values scale approximately as <r li ,a 1 oc ogf^ ' 7-0 ' 8 at a 
fixed source redshift, and increase with the cosmological constant because of the resulting increase in the 
path length (Bernardeau et al. 1997; Jain & Seljak 1997). In all models, the rms values obtained from our 
ray-tracing experiments are significantly larger than the predictions of the linear method. These rms values 
should be compared preferentially to the predictions of the filtered nonlinear method, because the effect of 
lensing by structures smaller scale than the softening length is not included in our experiments. The rms 
values obtained from our experiments are slightly smaller than those predicted by the filtered nonlinear 
method for the E-dS model, and slightly larger for the O and A models. It is important to recall that the 



Table 3. The minimum values of the fluctuation of the lensing convergence, 8k 



Redshift of source z s E-dS O A 



1.0 -0.1243 -0.04504 -0.06567 

2.0 -0.3285 -0.1274 -0.2071 

3.0 -0.5500 -0.2203 -0.3724 
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Fig. 4. — Root-mean-squares of the lensing convergence and shear of images as a function of redshifts of 
source. The filled circles and filled triangle are for the lensing convergence and shear of images evaluated 
from the results of ray-tracing experiments, respectively. The solid curves represent predictions of the power 
spectrum approach with the nonlinear power spectrum by Peacock & Dodds (1996), dashed curves are those 
obtained from a top-hatted density field with a smoothing scale of R s = 0.1/i _1 Mpc, and long-dashed curves 
are the predictions of linear perturbation theory, (a) E-dS model; (b) O model; (c) A model. 

same trends were shown in the two-point correlation functions of the matter (Figure [l]). These two trends 
are consistent with each others at least qualitatively; the correlation function and the rms values are either 
both overestimated or both underestimated by the nonlinear predictions, depending upon the cosmological 
model. This suggests that the rms values obtained from our experiments are consistent at least qualitatively 
with the predictions of the filtered nonlinear method. 

These various comparisons show that rms values of the lensing convergence and shear are identical 
within the uncertainties inherent to the numerical method. They also show that the power spectrum 
approach with the nonlinear power spectrum provides good estimates for the rms values of the lensing 
convergences and shears. 



3.2. Lensing Magnification 

Next, we study the statistics of the lensing magnifications. In order to obtain meaningful statistics, 
and be able to study properties of light rays with a large magnification, which are very rare among random 
lines of sight, we perform 10 experiments for each model. 

In Figure ^, we plot the probability distributions of image magnifications of a point-like source. 
For each model, the top panel shows the differential distribution p(n), and the bottom panel shows the 
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Fig. 5. — Probability distributions of image magnifications. The solid lines are for the source redshift z s = 3, 
dashed lines are for z s = 2, and dotted lines are for z 8 = 1. (a) E-dS model; (b) O model; (c) A model. 
For each model, top panel shows the differential distribution and bottom panel shows the integral one 
P(> AO- 
integral distribution p(> /i). The probability distributions are normalized so as to satisfy the normalization 
constraint dfip([i) = 1. In the differential distributions, we find that the distributions have a strong 
peak located at values of fi slightly less than 1, and are significantly skewed toward large magnifications. 
In the case of the E-dS model, (Figure ||a), a power-law tail of the form oc /i~ 2 appears at high 
magnifications, for sources located at redshifts z s — 2 and 3. This power-law behavior is a general property 
of the lens equation for a point-like source (see, e.g., SEF, chapter 11). The appearance of this power-low 
tail indicates the existence of caustics. 

The magnification distributions shown in Figure |5| can be interpreted as the "transfer function" of the 
matter in the universe (Wambsganss et al. 1998). In other words, any intrinsic luminosity function of, say, 
quasars will be folded with this magnification distribution, and what we measure as the observed quasar 
luminosity function will be the convolution of the intrinsic luminosity function with this transfer function 
of the universe (Turner 1980; Avni 1981; Canizares 1982, Peacock 1982; Victri & Ostriker 1983; Vietri 
1985; Ostriker & Vietri 1986; Schneider 1992). Therefore, it is essential to study first the magnification 
distribution in order to compute the magnification bias accurately. We shall discuss this point in detail in 
§4. 
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Since gravitational lensing not only magnify sources but also causes a distortion of their area on the 
sky, the lines of sight used in our experiments, which are randomly distributed on the image plane are 
not randomly distributed on the source plane (Ehlers & Schneider 1986). The small area of a source on 
the source plane is enlarged by a factor of /i on the image plane. Therefore if one averages the lensing 
magnifications over images randomly distributed on the image plane, the resulting average magnification 
(/j.) will be larger than unity. Taking the effect of area distortion into account, the probability distribution 
of the lensing magnifications of sources, i.e., the probability distribution calculated from the magnifications 
of sources randomly distributed on a source plane, is given by 



Notice that, by definition, the mean magnification over sources automatically becomes unity, 



Using the probability distributions shown in Figure |5|, we calculated the dispersions of the magnification 
distributions, 



where the probability distribution functions p s {[i) are sampled in bins equally spaced in intervals of 0.01 
in log/Lt, in the range 0.1 < /i < 100. The results are shown as filled circles in Figure El The solid curves 
represent the predictions of the power spectrum approach, using the nonlinear power spectrum of Peacock 
& Dodds (1996), and the dashed curves represent predictions based on the top-hat filtered density field 
with a smoothing scale of R s = 0.1/i _1 Mpc. The dispersions in the E-dS model are larger than those in the 
O and A models by about a factor of 3. In all models, the magnification dispersion is larger than 0.1 for 
a source redshift z s = 3. The effect of magnification dispersion on observations of sources will clearly be 
important for such high redshifts, especially if the density parameter f2o is large. 

Figures 6b and 6c show that, for the O and A models, the dispersions obtained from the experiments 
are in good agreement with the predictions of the nonlinear method, but are larger that the predictions of 
the filtered nonlinear method, as in the cases of lensing convergences and shears. This is again consistent 
with the trend found in the two-point correlation function at least qualitatively; i.e., the fitting formula of 
the power spectrum of Peacock & Dodds (1996) underestimates the two-point correlation functions for the 
O and A models. Figure |]a shows that for the E-dS model, the rms values agree with the predictions of 
the filtered nonlinear method when the source redshift is smaller than 1. However, the growth of ct m with 
the source redshift is more rapid than the filtered nonlinear method predicts. This rapid growth reflects 
the appearance of the power-law tail in the magnification distribution. In the power spectrum approach, a 
perturbative treatment is adopted for calculating the lensing magnification, and thus high magnification 
events are not properly taken into account. Therefore, the dispersion of lensing magnification estimated 
from the power spectrum approach should be regarded as a lower limit. 

Notice that the integral in equation ( |2l| ) diverges once the power-law behavior at the high magnification 
appears. In our calculations, this divergence is prevented by the fact that we have a finite number of 
experiments: the probability p s eventually drops to zero at some value of n, simply because no experiments 
have produced larger values. In the real universe, the integral does not diverge because the power-law tail 
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Fig. 6. — The dispersions in the lensing magnification of sources <7 M versus source redshift. The filled circles 
are evaluated from the results of the ray-tracing experiments. The solid curves represent predictions of the 
nonlinear method, dashed curves are those of the filtered nonlinear method, (a) E-dS model; (b) O model; 
(c) A model. 

is effectively cut off at large magnifications. This is caused by the following two effects. First, astrophysical 
sources are extended, and their magnifications (given by the surface brightness-weighted point-source 
magnification over the solid-angle area of the source) remain finite (e.g. Bontz 1979; Schneider 1987a; see 
also SEF, chapter 12). Second, even a point source can only be magnified by a finite value, because the 
geometrical optics approximation is inadequate near critical curves. A more accurate treatment, based on a 
wave-optics description of lensing, shows that the magnification is always finite (e.g., Ohanian 1974; 1983; 
Bliokh & Minakov 1975; Nakamura 1998; also see SEF, chapter 7). 

We should note here that, as was mentioned above, the rms values are sensitive to the normalization of 
the density contrast field, as, i.e. roughly a^ oc <t k ,ct 7 oc as (Bernardeau et al. 1997; Jain & Seljak 1997; 
Nakamura 1997). Thus the large rms values and the appearance of the power-law tail in magnification 
distribution in E-dS model are partly due to the high normalization in the model. Observations of the local 
cluster abundance also place constraints on the value of as which suggest low normalizations a 8 ~ 0.5 to 
0.6 for the fiat model with Ct = 1 (e.g., Eke, Cole & Frcnk 1996; Kitayama & Suto 1996, 1997; Viana & 
Liddlc 1996). The E-dS model might be, therefore, considered as an extreme case. 
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4. EFFECT OF THE MAGNIFICATION BIAS ON QUASAR LUMINOSITY 

FUNCTIONS 

Strictly speaking, almost all cosmological observations are under the influence of the lensing 
magnification, but there are very large differences in degree; the majority of field distant sources (say, 
redshift z s > 1) might be demagnified slightly (Wambsganss et al. 1998), whereas strongly lensed sources 
such as multiply- imaged quasars and giant luminous arcs, which are relatively very rare, are highly 
magnified (see e.g. SEF; and for a recent review, Mellier 1999). 

In this section, we study the effect of the magnification bias on quasar luminosity functions, using the 
magnification distributions obtained from the ray-tracing experiments described in §§2 and 3. Although 
magnification bias is not a phenomenon restricted to quasars alone, we focus our attention on these sources, 
as previous authors have done (e.g., Schneider 1992), because lensing effects are most relevant to them, due 
to their inferred compactnesses and large cosmological distances. 

We use the conventional definition of the luminosity function <f>(L,z) of quasars, as being the comoving 
number density of quasars at a redshift z with a luminosity L. We designate by <I?i n t(L,.z) the intrinsic 
luminosity function of quasars, which is the actual luminosity function, and by $obs(^,^) the observed 
luminosity function, which is the luminosity function inferred from observations. Observers compute 
<&obs{L, z) by counting sources in redshift bins, and estimating their luminosities from their apparent 
brightnesses and the luminosity distance Dl(z) corresponding to each redshift bin. In general, <J? bs(£, 2) 
and $i n t(L, z) differ in presence of gravitational lensing. A source believed to have a luminosity L might 
actually be a source with luminosity L/p, which is magnified]^] by a factor p. If we neglect any other 
uncertainty in the determination of ^obs, then <& b s (L, z) and (&i n t(L, z) are related by 

$ ohs {L,z)=J^ dfx^ int f^,z\p s (p,z) (22) 

(SEF, chapter 12). In the absence of lensing, that is, in a smooth Friedmann universe, p s (p>, z) is equal to a 
Dirac 5-function 5d(p< — 1), and equation ( p2[ ) reduces to ^obs = ^int- We define the bias parameter 6q as 

_ & obB (L,z) 

<Pint(£, Z) 

To compute the magnification bias, we need two ingredients: the magnification probability p s (p, z) and the 
intrinsic luminosity function $; nt (L,z). The former is provided by our ray-tracing experiments. For the 
latter, we consider two empirical models: a single power-law model, and a double power-law model. 



4.1. Single Power-law Model 

Let us first consider an intrinsic luminosity function described by a single power-law model, 

—ct 



<!>;„, (I.:.) <I> : (:)(-^) . (24) 



5 or demagnified, since fi can be either larger or smaller than unity. 
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where the amplitude <&*(z) depends on redshift but not luminosity. The power index a of the luminosity 
function is uncertain, but is believed to be in the range a ~ 1.2 — 1.6 at the faint end of the luminosity 
function, and a ~ 3 — 4 at the bright end (Boyle et al. 1988; La Franca & Cristiani 1997). Hence, a single 
power-law model is not a particularly good approximation for & ni t- However, this model is interesting from 
a theoretical viewpoint. 

We substitute equation ( |24| ) in equation (^||), and get 



Notice that &q is a function of z, but not L. Hence, in the case of the single power-law model, the 
logarithmic slope of the luminosity function is unchanged by the magnification bias, but the amplitude is 
multiplied by a redshift-dependent factor, and can cither increase or decrease, depending on the power 
index a and the probability distribution p s (ji,z). Notice also that in the particular case of a power-law 
index a = 1, we get $ bs = 'I'int- The observed luminosity functions in a smooth Friedmann universe and in 
an inhomogeneous Friedmann universe are then equal, in spite of the presence of lensing magnification. 

We computed the bias parameter &q as a function of the power-law index a, using the probability 
distribution functions derived from our ray-tracing experiments. As in §3.2, we limited the integration 
range to 0.1 < fi < 100. We considered three particular source redshifts, z s = 1, z s = 2, and z s — 3. The 
results are plotted in Figure 0. The bias strongly depends on the power index a. Since the actual quasar 
luminosity function becomes steeper as the luminosity increases, the effect of magnification bias should be 
more important at the bright end of the luminosity function than at the faint one. Figures [?]b and |?]c show 
that the biasing effect is moderate for the O and A models, and is at most 10% even for sources at z s = 3. 
On the other hand, the biasing is considerable for the E-dS model. Even for sources at moderate redshift, 
z s = 1, the bias parameter 6q is 4 for a — 4. For sources at z s = 2 and 3, the bias parameter strongly grows 
with a for the E-dS model. This is mainly due to the presence of the power-law tail in the magnification 
distributions that we saw in Figure ^|a. When the power-law tail appears in the magnification probability, 
the bias parameter scales roughly as 6q oc where /x max is the upper value of the integration over ^i, in 

our case /i ma x = 100. This implies that the magnification bias diverges in the limit of large magnifications, 
corresponding to low luminosities, when a > 2. This divergence is essentially an artifact of the single 
power-law model used here. The actual luminosity function flattens at the faint end, and the power-law 
index drops below 2, thus eliminating the divergence. Furthermore, the finite size of the source size becomes 
an important effect for very highly magnified sources (say \i > 10 3 , Schneider 1987c, 1992). Schneider & 
Weiss (1988b) and Schneider (1992) have pointed out that the magnification distribution is effectively cut 
off like p(/i) oc (At c //i) 6 for sources with a finite extend, where /j, c is a function of the source extend and the 
mass of a lens. In the case we considered (i.e. distant quasars lensed by large-scale structure), [i c is larger 
than 10 3 (see, e.g., SEF, chapter 12). We do not need to take this effect into account here, because the 
upper limit of the integration, /i max , is smaller than \i c . 




(25) 



In this case, the bias parameter is 




(26) 
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Fig. 7. — The bias parameter &q as a function of a 
power index a. The intrinsic luminosity function is 
taken by the single power-law form, equation (]24|). 
The solid lines are for the source redshift z s = 1, 
dashed lines are for z s = 2, and long-dashed lines 
are for z s = 3. (a) E-dS model; (b) O model; (c) A 
model. 
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Fig. 8. — The bias parameter 6q as a function of 
a luminosity normalized by its characteristic value, 
L/L*. The intrinsic luminosity function is taken 



by the double power-law form, equation (27), with 
ax = 1.5, and c*2 = 3.5. The meanings of lines are 
the same as in Figure (a) E-dS model; (b) O 
model; (c) A model. 
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4.2. Double Power-law Model 

Since the power-law index of the luminosity functions appears to take different values at the bright 
and faint ends, a single power law model, such as the one considered in §4.1, does not provide a very good 
fit, for any value of the power-law index a. A much better fit can be obtained by considering a double 
power-law model of the form 

®*( z ) 

$ int(£, z) - (L/L , )Q1 +(jL/L *)"2> (27) 

where ol\ and «2 are the logarithmic slope of the faint and bright end sides of the luminosity function, 
respectively. We choose ot\ = 1.5, and ot^ — 3.5 which are typical values (Boyle et al. 1988; La Franca & 
Cristiani 1997). 

In Figure |^, we plot the biasing parameter 6q in units of its characteristic value L* . Figures |§|b and 
||c show that the magnification bias is moderate for O and A models. The bias enhances the luminosity 
functions by about 5% at most at the bright end, and has little effect (below 1%) at the faint end. On the 
other hand, the effect is significant in the E-dS model, especially for sources at high redshift, z s > 2. At the 
present time, high-z quasars seems to be homogeneously sampled to 2 or 3 magnitude brighter than their 
characteristic value M* (which corresponds to a luminosity L ~ 10L*, e.g., La Franca & Cristiani 1997). 
Figure ||a shows that, for the E-dS model, the bright end of luminosity function (at L ~ 10L*) is enhanced 
by a factor larger than 2 for such high-z quasars. In this case, the observed luminosity function of bright 
quasars is strongly biased by the lensing magnification effect. 

As an illustrative example, we plot in Figure ^, for the E-dS model, the observed luminosity function 
as a function of the luminosity in units of its characteristic value L* , together with the intrinsic luminosity 
function, for sources at redshifts z s = 2 (dashed line) and at z s — 3 (long-dashed line). The lensing 
bias flatten the slope of the luminosity functions at the bright end of the luminosity function, above the 
characteristic luminosity L* . If wc fit the observed luminosity functions using the same double power-law 
model (eq. |27])]) as for the intrinsic luminosity function, we find that the values of L» and u\ are the 
same, but the "observed" power-law index ai takes the values ck2,obs = 3.0, and a2,obs = 2.8 for z s = 2 and 
3, respectively, compared to 012 — 3.5 for the intrinsic luminosity function. Since the observed power-law 
index depends on the source redshift, it is essential to take into account the lensing magnification bias 
when studying the evolution of the intrinsic quasar luminosity functions, in order to untangled the effects 
of magnification bias and genuine luminosity and/or number evolution. 



5. SUMMARY AND CONCLUSION 

We have studied the statistical properties of gravitational lensing by large-scale structure, for three 
different COBi^normalizcd cosmological models. We used a P 3 M TV-body code to simulate the formation 
and evolution of large-scale structure, and then used the multiple lens-plane algorithm to follow light rays 
propagating through the inhomogeneous matter distribution in the universe. For each model, we followed 
the propagation of 1.1 x 10 7 light rays, from the observer up to a source redshift of 3, and computed the 
evolution of the Jacobian matrix along each ray. Having such a large number of light rays enabled us to 
compute, as a function of source redshift, the distributions and rms values of the lensing convergences, 
shears, and magnifications, to high accuracy. We compared these results with various analytical estimates 
based on the power spectrum method. Finally, by combining the magnification probability obtained from 
the experiments with empirical models for the intrinsic luminosity function, of quasars, we computed the 
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Fig. 9. — The observed luminosity function at the redshifts z = 2 (the dashed line), and z = 3 (the 
long-dashed line) are plotted as a function of the luminosity normalized by its characteristic value, L/L*, 
together with the intrinsic luminosity function (the solid line) . The intrinsic luminosity function is taken by 
the double power-law form, equation ( |27| ) with the power indices ct\ — 1.5, and Q2 = 3.5. The background 
cosmology is taken by E-dS model. 

magnification bias. Our main results can be summarized as follows: 

1. The rms values of the lensing convergences and shears, which are expected to be equal according 
to analytical arguments, were found to be nearly equal, the difference being attributed to the finite 
resolution of the numerical methods used for the simulations. These values were consistent with 
analytical predictions based on the power spectrum approach with nonlinear evolution of the power 
spectra. 

2. The dispersion of the lensing magnifications strongly depends on the amplitude of the density 
fluctuations of the matter and thus depends on the density parameter f2 - The dispersion becomes 
larger than 0.1 at z s = 1 for the E-dS model, and at z s — 3 for the O and A models. The magnification 
distributions are considerably skewed toward high magnification. In particular, a power-law tail 
appears at large magnifications in the E-dS model, for sources at redshifts z s > 1. 

3. We have compared the statistics of the lensing magnification with the nonlinear predictions of the 
power spectrum approach. We found that the power spectrum approach with the nonlinear power 
spectrum correctly predicts the value of the dispersion of the magnification in the absence of the 
power-law tail. However, once the power-law tail appears, the predicted values become inaccurate, 
because the dispersions grow more rapidly with the source redshift than the nonlinear predictions. 
These predictions, therefore, should be regarded as a lower limit. 
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4. Wc studied the magnification bias on quasar luminosity functions. We found that the lensing 
magnification bias strongly depends on the slope of the luminosity function and on the density 
parameter Qq. If the mean matter density is as high as in the E-dS model, the biasing effects can 
be significant, especially at the bright end of the luminosity function where its slope is very steep. 
Moreover, since quasar luminosity functions becomes steeper with the luminosity, the lensing bias will 
flatten the effective slope of an observed luminosity functions at the bright side. 

The power-law tail did not appear in the O and A models because the large-scale structure in these 
models is not as evolved as in the E-dS model. However, the existence of strongly magnified images, such 
as multiply imaged quasars or galaxies and giant luminous arcs, suggests that whatever the background 
cosmological model is, there is a power-tail in the actual magnification distribution. Therefore, in order to 
study the effect of lensing magnification bias on the observed luminosity function, one must carefully take 
into account the influence of this power-law tail. It is also important to emphasize that, as we have stated 
in §3.2, the large lensing effects found in the E-dS model are partly due to the high normalization used for 
that model. Therefore, the E-dS model should be considered as an extreme case. 

We have investigated the magnification bias resulting from the gravitational lensing in universes with 
realistic matter distributions. In order to study quantitatively the biasing effects on the observed quasar 
luminosity functions, one has to take into account additional selection effects such like the finite extend of a 
source, which have not been considered in this paper. Such effects will be studied in future works. 
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Grants- in- Aid by the Ministry of Education, Science, Sports and Culture of Japan (09640332). HM was 
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In this Appendix, we present the analytic formulae for calculating the variances of the lensing 
convergences, shears, and magnifications. Those formulae have been derived by Kaiser (1992, 1998), 
Villumsen (1996), Bernardeau et al. (1997), Frieman (1997), Jain & Scljak (1997), and Nakamura (1997). 
We refer the reader to these references for details, and only present the final expressions. 

In the weak gravitational field limit, the variance of the lensing convergence and shear can be estimated 
using both the Born approximation (Schneider et al. 1998) and Limber's equation in Fourier space (Kaiser 
1992; 1998). These variance, which turn out to be equal, are related to the density power spectrum, by 



where x is the comoving distance, -D(xi,X2) is the standard angular diameter distance between \i an d Xi-, 
and a is the Robertson- Walker scale factor normalized to be unity at the present. The relationship between 
the comoving distance and the redshift z can be derived from the Friedmann equation (Jain & Seljak 1997): 



A. 



CALCULATION OF THE VARIANCES 
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The function I(v) appearing in equation (Al) in an integral over the power spectrum. It is the particular 
form of that integral that distinguishes the various approximations. In linear theory, this function is given 

by 

/■OO 

I(v) = / dkkPL[k,a(v)], (A3) 
Jo 

where Pl is the linear power spectrum, which is a function of the wavenumber k and the scale factor a, 
and the dependence of the comoving distance v enters only through the scale factor. In the nonlinear 
approximation, we use instead 

/•oo 

I(v) = / dkkP Nlj [k,a{v)}, ( A4 ) 
Jo 

where Pnl is the nonlinear power spectrum, for which Peacock & Dodds (1996) provide a fitting formula. 
Finally, for the top-hat filtered density field approximation, we use 

/>oo 

I(v,R s )= dkkP NIj [k,a(v)}W% H (kR s ), (A5) 
Jo 

where Wth(x) = 3(sinx — xcosx)/x 3 is the Fourier transform of the top-hat window function with 
smoothing scale R s . The expression of the variance of the lensing magnifications have been derived by 
Frieman (1997) and Nakamura (1997), who found that er^ = 4^2 
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